/*
 * Basic component library.
 *
 * Copyright (c) 1993-2010 Leo Spiekman, Xaveer Leijtens
 *
 * This file is part of libMask.
 *
 * libMask is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * libMask is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with libMask.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "mask.h"

#define VERYSMALL   (1.0e-12)

/* Maximum number of point in a polyline or polygon */
#define MAXPTS	    (200)

extern void npolyline(double*, double*, int, double);
extern void npolygon(double*, double*, int, int, int);

/* Draws a curved section of two waveguides with a conductor width of w and
 * a gap of g. For a curve to the right, both a and radius r should be
 * positive, for a turn to the left, a and r should be negative
 */
void cw2(double r, double a, double w, double g)
{
  push_cp();
  offset((g + w) / 2);
  cw(r - (g + w) / 2, a, w);
  pop_cp();
  offset(-(g + w) / 2);
  cw(r + (g + w) / 2, a, w);
  offset((g + w) / 2);
}

void curve_rect(double w, double r, double a, int n_wg)  {

/* Defintions */

double l;                           /* Length of one block */
double dl;                          /* Additional length */
double an;                          /* single angle of n-hook*/

int i;                              /* counter for the tapers */

  an=RAD(a/(n_wg-1));
  l = fabs(2*r*tan(an/2));
  dl = tan(an/2)*w/2;  
  if (a<0)
    dl = -dl;
  rect(l/2+dl,w,4,6);
  move(-dl,0,0,NOFLIP);
  move(0,0,an,NOFLIP);
  move(-dl,0,0,NOFLIP);
  for (i=0;i<n_wg-2;i++) {
    rect(l+2*dl,w,4,6);
    move(-dl,0,0,NOFLIP);
    move(0,0,an,NOFLIP);
    move(-dl,0,0,NOFLIP);
  }
  rect(l/2+dl,w,4,6);
}

void cwg_curve_rect(double w, double g, double r, double a, int n_wg_i, int n_wg_o)  {

/* Defintions */

double dx;                           /* Length of the curve */
double dy;                          /* y position of the center of the curve at the end*/
double w_tot;                       /* total width of the block*/

  w_tot=w+2*g;
  dx=-r*sin(RAD(a));
  dx=(a<0)?dx:-dx;
  dy=-r*(1-cos(RAD(a)));
  dy=(a<0)?dy:-dy;
  movept(dx,w_tot-g,4,7);
  push_cp();
    if (a<0)
      curve_rect(g,r+(w+g)/2,a,n_wg_o);
    else 
      curve_rect(g,r-(w+g)/2,a,n_wg_i);
  pop_cp();
  movept(dx,w_tot-g,7,1);
  push_cp();
    if (a<0)
      curve_rect(g,r-(w+g)/2,a,n_wg_i);
    else
      curve_rect(g,r+(w+g)/2,a,n_wg_o);
  pop_cp();
  movept(dx,w_tot-g,1,6);
  move(0,dy,0,NOFLIP);
  move(0,0,RAD(a),NOFLIP);

}

/* This implementation does not work with ADS and MDS. See comment in egs.c
 * for egsarc2().
 */
void cw_35ph(double r, double a, double w)
{
  double x,y;
  double a1,a2;
  double s=M_PI/2;

  r = copysign(r,a); /* r gets the sign of a */
  s = copysign(s,a);

  /* center point */
  offset(r);
  x=cp.x;
  y=cp.y;

  /* first angle */
  push_cp();
  rotate (-s);
  a1 = cp.a;
  pop_cp();

  rotate(a);

  /* second angle */
  push_cp();
  rotate (-s);
  a2 = cp.a;
  pop_cp();

  if ((fabs(a1 - a2)-fabs(a)) > M_PI/2) {
    if (a1 <= 0) a1+=M_PI*2;
    if (a2 <= 0) a2+=M_PI*2;
  }

  offset(-r);

  if (a1 > a2) mskarc2(clayer,w,fabs(r),a2,a1,x,y);
  else mskarc2(clayer,w,fabs(r),a1,a2,x,y);

}

/* this implementation does not work with ADS and MDS. See comment in egs.c
 * for mskarc2().
 * This routine has NOT been thoroughly tested.
 */
void cw_35ph_b(double r, double a, double w)
{
  double x,y;
  double a1,a2;
  double s=M_PI/2;

  r = copysign(r,a); /* r gets the sign of a */
  s = copysign(s,a);

  /* center point */
  offset(r);
  x=cp.x;
  y=cp.y;

  a1 = map_0_2pi(cp.a + M_PI/2);
  rotate(a);
  a2 = map_0_2pi(cp.a + M_PI/2);

  offset(-r);

  if (a1 > a2) mskarc2(clayer,w,fabs(r),a2,a1,x,y);
  else mskarc2(clayer,w,fabs(r),a1,a2,x,y);
}

/* Make curve with polylines (by not setting mode CWARC).
 * Approximate the arc by sections of straights connecting points on an arc
 * with radius Re = R / cos(da). The first and last steps are at da/2, the
 * rest are da apart.
 * The maximum deviation has to be less than for the polygon version,
 * because that routine uses an adapted radius that prevents the systematic
 * error present in this approach.
 */
void cw_polyline(double r, double a, double w)
{
  const double eps = getlayer_accuracy(currentlayer());
  double astart, da, ang;
  double Re, R; /* Effective radius, fabs(r) */
  int i, k, n, nstep, nsection;
  double x[MAXPTS], y[MAXPTS];

  R = fabs(r)+w/2; /* Outside curve radius */
  r = copysign(r, a);
  astart = copysign(M_PI/2, a);
  /* Number of steps needed for the required accuracy, calculated from the
   * outside of the curve. This determines the number of sections.
   */
  da = 2 * acos(R/(R + eps));
  nstep = (int)(fabs(a) / da) + 2; /* 2: 1 for rounding, 1 for endpoint */
  /* A total of MAXPTS (200) is allowed by GDSII, but we split already
   * here.
   */
  nsection = nstep / (MAXPTS-2) + 1; /* The number of sections */
  for (k=0; k<nsection; k++) {
    da = 2 * acos(R/(R + eps));
    n = (int)(fabs(a/nsection) / da) + 2;
    da = a / nsection / (n-1);
    Re = fabs(r) / cos(da/2);
    x[0] = 0; /* first point */
    y[0] = 0;
    for (i=1; i<n; i++) {
      ang = astart + da/2 - i*da;
      x[i] = Re * cos(ang);
      y[i] = -r + Re * sin(ang);
    }
    ang = astart - a/nsection; /* last point */
    x[n] = fabs(r) * cos(ang);
    y[n] = -r + fabs(r) * sin(ang);
    npolyline(x, y, n+1, w);
    rotate(a/nsection/2);
    skip(2*r*sin(a/nsection/2));
    rotate(a/nsection/2);
  }
}

/* Make curve with polygons (by setting mode POLYGON).
 * Approximate the arc by sections of straights connecting points on an arc
 * with radius Re = R * da / sin(da). The first and last steps are at da/2,
 * the rest are da apart. The effective radius chosen like this, prevents
 * systematic errors on the length of the curve.
 */
void cw_polygon(double r, double a, double w)
{
  const double eps = getlayer_accuracy(currentlayer());
  double astart, da, ang;
  double Re, Ro, Ri; /* Effective, outer/inner radii */
  int i, k, nstep, nsection, no, ni;
  double x[MAXPTS-1], y[MAXPTS-1], xc, yc;

  Ro = fabs(r) + w/2; /* Outside curve radius */
  Ri = fabs(r) - w/2; /* Inside curve radius  */
  if (Ri < 0) {
    fprintf(stderr, "Waveguide half-width (%.3f um) larger than radius (%.3f um)\n", w*0.5e6, r*1e6);
    Ri = 0;
  }
  r = copysign(r, a);
  astart = copysign(M_PI/2, a);
  xc = 0;
  yc = -r;
  /* Number of steps needed for the required accuracy, calculated from the
   * outside of the curve. This determines the number of sections.
   */
  da = 2 * acos(Ro/(Ro + eps));
  nstep = (int)(fabs(a) / da) + 2; /* 2: 1 for rounding, 1 for endpoint */
  /* A total of MAXPTS (200) is allowed by GDSII, but we split already
   * here.
   */
  nsection = nstep / ((MAXPTS-1)/2-2) + 1; /* The number of sections */
  for (k=0; k<nsection; k++) {
    /* Outside curve */
    da = 2 * acos(Ro/(Ro + eps));
    no = (int)(fabs(a/nsection) / da) + 2;
    da = a / nsection / (no-1);
    Re = Ro * da/2 / sin(da/2);
    x[0] = xc; /* first point */
    y[0] = copysign(w/2, a);
    for (i=1; i<no; i++) {
      ang = astart + da/2 - i*da;
      x[i] = xc + Re * cos(ang);
      y[i] = yc + Re * sin(ang);
    }
    ang = astart - a/nsection; /* last point */
    x[no] = xc + Ro * cos(ang);
    y[no] = yc + Ro * sin(ang);
    /* Inside curve */
    da = 2 * acos(Ri/(Ri + eps));
    ni = (int)(fabs(a/nsection) / da) + 2;
    da = a / nsection / (ni-1);
    Re = Ri * da/2 / sin(da/2);
    x[no+1] = xc + Ri * cos(ang); /* first point backwards */
    y[no+1] = yc + Ri * sin(ang);
    for (i=1; i<ni; i++) {
      ang = astart + da/2 - (ni-i)*da;
      x[no+1+i] = xc + Re * cos(ang);
      y[no+1+i] = yc + Re * sin(ang);
    }
    x[no+1+ni] = xc;
    y[no+1+ni] = -copysign(w/2, a);
    npolygon(x, y, no+1, no+1, no+ni+2);
    rotate(a/nsection/2);
    skip(2*r*sin(a/nsection/2));
    rotate(a/nsection/2);
  }
}

/* This routine uses three points to define the arc */
void cw_3pts(double r, double a, double w)
{
  double x1, x2, x3, y1, y2, y3;

  r = copysign(r, a);		/* r gets the sign of a */

  x1 = cp.x;
  y1 = cp.y;
  rotate(a / 4);
  skip(2 * r * sin(a / 4));
  x2 = cp.x;
  y2 = cp.y;
  rotate(a / 2);
  skip(2 * r * sin(a / 4));
  x3 = cp.x;
  y3 = cp.y;
  rotate(a / 4);

  mskarc(clayer, w, x1, y1, x2, y2, x3, y3);
}
void cw(double r, double a, double w)
{
  double a_extra;
  int state;

  if (fabs(a) < VERYSMALL || fabs(r) < VERYSMALL) return;
  /* overlap */
  state = pen_up();
  if (state == PENUP) {
    /* pen was already up, nothing to write, just move the pointer */
    cw_polyline(r, a, w);
    return;
  }
  /* pen was down */
  a_extra = copysign(cw_overlap / r, a);
  if (fabs(a_extra) > VERYSMALL) {
    rotate(RAD(180));
    cw_polyline(r, -a_extra, w); /* back off a little */
    rotate(RAD(180));
  } else {
    a_extra = 0;
  }
  pen_down();
  if (mode & CWARC) {
    if (mode & CWANGLE)
      cw_35ph(r, a + 2 * a_extra, w);
    else
      cw_3pts(r, a + 2 * a_extra, w);
  } else if (mode & POLYGON) {
    cw_polygon(r, a + 2 * a_extra, w);
  } else {
    cw_polyline(r, a + 2 * a_extra, w);
  }
  /* move the pointer to the proper position */
  pen_up();
  if (fabs(a_extra) > VERYSMALL) {
    rotate(RAD(180));
    cw_polyline(r, -a_extra, w);	/* Go back to proper position */
    rotate(RAD(180));
  }
  pen_down(); /* restore pen state */
}

void corner90(double a, double w)
{
  double x1,y1,x2,y2;

  w=(a>0)?w:-w;
  offset(-w/2);
  x1=cp.x;
  y1=cp.y;
  offset(w);
  x2=cp.x;
  y2=cp.y;
  rotate(a);
  offset(-w);
  mskpolygon3(clayer,x1,y1,x2,y2,cp.x,cp.y);
  offset(w/2);
}

void circle(double w, double r)
{
  double x1,x2,y1,y2;

  if (mode & CWARC) {
    x1=cp.x;
    y1=cp.y;
    x2=x1+r;
    y2=y1;
    mskcircle(clayer,w,x1,y1,x2,y2);
  } else if (mode & POLYGON) {
    push_cp();
    offset(r);
    cw_polygon(r, -M_PI, w);
    cw_polygon(r, -M_PI, w);
    pop_cp();
  } else {
    push_cp();
    offset(r);
    cw_polyline(r, -M_PI, w);
    cw_polyline(r, -M_PI, w);
    pop_cp();
  }
}

/* filled circle */
void disc(double r)
{
  const double eps = getlayer_accuracy(currentlayer());
  int i, j, nstep, nsection;
  double da, xi, yi, xc, yc, step;

  r = fabs(r);
  /* Number of steps needed for the required accuracy. */
  da = 2 * acos(r/(r + eps));
  nstep = (int)(2 * M_PI / da) + 1; /* for rounding */

  /* A total of MAXPTS (200) is allowed by GDSII, but we split already
   * here.
   */
  if (nstep > MAXPTS-2) { /* more than one section */
    nsection = nstep / (MAXPTS-2) + 1; /* number of sections */
    nstep = nstep / nsection + 1; /* per section */
  } else {
    nsection = 1;
  }
  da = 2 * M_PI / (nstep * nsection);
/*  r = fabs(r) / cos(da/2); */
  r = r * da/2 / sin(da/2);

  push_cp();
  step = 2 * r * sin(da/2);
  xc = cp.x; /* center coordinates */
  yc = cp.y;
  move(0, r, -da/2, 0);
  for (j=0; j<nsection; j++) {
    xi = cp.x;
    yi = cp.y;
    mskpolygonstart(clayer);
    if (j > 0) mskpolypoint(clayer,xi,yi);
    for (i=0; i<nstep; i++) {
      mskpolypoint(clayer,cp.x,cp.y);
      move(step, 0, -da, 0);
    }
    mskpolypoint(clayer,cp.x,cp.y);
    if (nsection > 2) /* include center in polygon (piece of cake) */
      mskpolypoint(clayer,xc,yc);
    if (nsection > 1) /* close path: go back to start point */
      mskpolypoint(clayer,xi,yi);
    mskpolygonend(clayer);
  }
  pop_cp();
}

double adiabend(double Rmini, double theta, double width, int N)
{
  double dTheta = theta / (2 * N);
  double L = 0, Ri;
  int i;

  for (i = 0; i < N; i++) {
    Ri = Rmini + (Rmini * 7 * pow((N - 1 - i), 7)) / pow((N - 1), 7);
    cw(Ri, dTheta, width);
    L = L + dTheta * Ri;
  }
  for (i = 0; i < N; i++) {
    Ri = Rmini + (Rmini * 7 * pow(i, 7)) / pow((N - 1), 7);
    cw(Ri, dTheta, width);
  }
  return (2 * L);
}
